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ABSTRACT 

Performing one-dimensional hydrodynamical calculations coupled with nonequilibrium 
processes for hydrogen molecule formation, we pursue the thermal and dynamical evo- 
lution of filamentary primordial gas clouds and attempt to make an estimate on the 
mass of population III stars. The cloud evolution is computed from the central proton 
density n c ~ 10 2-4 cm -3 up to ~ 10 13 cm~ 3 . It is found that, almost independent of 
initial conditions, a filamentary cloud continues to collapse nearly isothermally due to 
H2 cooling until the cloud becomes optically thick against the H2 lines (n c ~ 10 10 ~ n 
cm -3 ). During the collapse the cloud structure separates into two parts, i.e., a denser 
spindle and a diffuse envelope. The spindle contracts quasi-statically, and thus the line 
mass of the spindle keeps a characteristic value determined solely by the temperature 
(~ 800 K), which is ~ 1 x 10 3 M© pc" 1 during the contraction from n c ~ 10 5 cm~ 3 to 
10 13 cm -3 . Applying a linear theory, we find that the spindle is unstable against frag- 
mentation during the collapse. The wavelength of the fastest growing perturbation (A m ) 
lessens as the collapse proceeds. Consequently, successive fragmentation could occur. 
When the central density exceeds n c ~ 10 10 ~ n cm -3 , the successive fragmentation may 
cease since the cloud becomes opaque against the H2 lines and the collapse decelerates 
appreciably. Resultingly, the minimum value of A m is estimated to be ~ 2 x 10~ 3 pc. 
The mass of the first star is then expected to be typically ~ 3M , which may grow 
up to ~ 16M Q by accreting the diffuse envelope. Thus, the first-generation stars are 
anticipated to be massive but not supermassive. 



Subject headings: cosmology — galaxies: formation — hydrodynamics — ISM: clouds 
— stars: formation 



1. Introduction 



Based on the standard Big Bang nucleosynthesis, the first generation of stars should form from 
materials deficient in heavy elements. In the present-day galaxies, heavy elements or dust grains 
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provide the most efficient cooling mechanism, while the cooling process in primordial gas is likely to 
be governed by hydrogen molecules. Many authors hitherto have considered the formation processes 
of primordial stars from metal deficient gas (e.g., Matsuda, Sato, & Takeda 1969; Yoneyama 1972; 
Hutchins 1976; Silk 1977a, 1977b; Yoshii & Sabano 1979; Carlberg 1981; Struck-Marcell 1982a, 
1982b; Lepp & Shull 1983, 1984; Silk 1983; Palla, Salpeter, & Stahler 1983; Yoshii & Saio 1986; 
Shapiro k Kang 1987; de Araiijo & Opher 1989; Uehara et al. 1996; Haiman, Thoul, & Loeb 1996; 
Omukai et al. 1998). Such first-generation stars, say Population III, could play an important role 
in the early evolution of galaxies (e.g., Tegmark, Silk, & Blanchard 1994; Ostriker & Gnedin 1996) 
or the formation of massive black holes (Umemura et al. 1993). Also, they may be responsible for 
the chemical pollution of intergalactic medium which is recently inferred from metallic absorption 
in Lya forest seen in quasar light (Cowie et al. 1995; Songaila & Cowie 1996). It is thus important 
to study the thermal and dynamical evolution of the primordial gas clouds from which the first- 
generation stars would form. 

If the hydrogen gas were to remain purely atomic, the primordial gas would be cooled down to 
10 4 K due to the Lyman a lines. It is, however, difficult for the gas temperature to become lower, 
because hydrogen atoms are poor radiator in the lower temperature gas. Therefore, the Jeans mass 
of such a cloud, which is often referred to the characteristic mass in the star formation theory, 
becomes much greater than a typical stellar mass. In practice, hydrogen molecules provide key 
cooling mechanisms. In contrast to the molecule formation on dust grains in metal-rich interstellar 
gas, the formation of primordial hydrogen molecules can proceed through the gas phase reaction 
(Saslaw & Zipoy 1967; Peebles &Dicke 1968), 

e + H -» H - + hu (1) 

IT + H -» H 2 + e (2) 

and 

H+ + H -» FP+ + hu (3) 

H+ + H H 2 + H+. (4) 

At high density of n ^ 10 8 cm~ 3 , the molecules can also form through the three-body reactions 
(Palla et al. 1983), 

3H -> H 2 + H (5) 

and 

2H + H 2 -» 2H 2 . (6) 

Even if a relatively small fraction (~ 10~ 3 ) of the molecules form, they make a significant con- 
tribution to cooling through the rotational and vibrational transitions, so that the temperature 
of primordial gas can be reduced to a lower temperature than 10 4 K (e.g., Matsuda et al. 1969; 
Yoneyama 1972; Hutchins 1976; Palla, et al. 1983). Resultantly, the Jeans mass could decrease to 
stellar mass. However, although a lot of elaborate analyses have been made so far, the masses of 
Population III stars have not been well converged. Carlberg (1981) and Palla et al. (1983) have 
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shown the Jeans mass can go down to a level of <^ O.IMq. Yoshii & Saio (1986) derive the initial 
mass function (IMF) based upon the opacity-limited fragmentation theory (Silk 1977a, 1977b) and 
find the peak of the IMF around 4 — WM & . Uehara et al. (1996) have shown that the mini- 
mum masses of the first-generation stars are basically determined by the Chandrasekhar mass, say 
~ IMq (see also Rees 1976). In this paper, we reanalyze the formation of Population III stars by 
computing the collapse of filamentary clouds coupled with R~2 formation. 

In the bottom-up scenarios like a cold dark matter (CDM) model, it is expected that the 
overdense regions with the masses of 10 5 — 1O 7 M0 would first collapse at the redshift range of 
10 ^ z ^ 100 an d the first generation of stars would form there. The importance of the H2 cooling 
in the collapse of such primordial clouds has been stressed by several authors (e.g., de Araiijo & 
Opher 1989; Haiman et al. 1996; Susa et al. 1996; Tegmark et al. 1997). It has been found 
in these studies that the mass fraction of H2 can reach 10 -4 up to 10 -3 and the temperature of 
the cloud can be reduced to 10 2 — 10 3 K. However, most of the studies are restricted to highly 
simplified models such as homogeneous pressure-less and/or spherical collapses. In practice, the 
cloud contraction proceeds inhomogeneously. Since the cloud is more or less nonspherical, the 
deviation from spherical symmetry grows in time due to its self-gravity until a shocked pancake 
forms (e.g., Umemura 1993). Recently, several authors have studied the collapse of pregalactic 
clouds with multi-dimensional hydrodynamical simulations (e.g., Anninos & Norman 1996; Ostriker 
& Gnedin 1996). Unfortunately they still do not have enough spatial resolution as well as mass 
resolution to explore the star formation. 

The pancake could be gravitationally unstable against fragmentation (e.g., Umemura 1993; 
Anninos &: Norman 1996). It tends to fragment into thin filamentary clouds rather than spherical 
ones (Miyama, Narita, & Hayashi 1987a, 1987b; Uehara et al. 1996). The filamentary cloud can also 
fragment into smaller and denser cores (e.g., Larson 1985), in which consequently stars can form. In 
this paper we thus investigate the thermal and dynamical evolution of the filamentary primordial 
gas clouds. With using one-dimensional axisymmetric hydrodynamical scheme, we pursue the 
evolution in the range of more than eighth order of magnitude in density contrast to properly 
estimate the mass scale of the first-generation stars. 

In we describe our model clouds and the computational methods. Numerical results are 
presented in §|3]. In §|| we consider the fragmentation of a collapsing filamentary cloud and estimate 
the mass scale of the fragments and thereby the mass of the first-generation stars. §|5] is devoted to 
the conclusions. 



2. MODEL 
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2.1. Basic Equations 

To pursue the thermal and dynamical evolution of a filamentary primordial cloud, we employ 
a one-dimensional hydrodynamical scheme. We assume that the system is axisymmetric and that 
the medium consists of ideal gas. The adiabatic index, 7, is taken to be 5/3 for monatomic gas and 
7/5 for diatomic gas. We deal with the following nine species: e, H, H + , H~, H2, H;j~, He, He + , 
and He ++ , The abundance of helium atoms is taken to be 10 % of that of hydrogen by number. 

The basic equations are then given in the cylindrical coordinates (r, ip, z) by 

dp 1 d , 

d , . 1 9 , 9. dP dib , . 

+ + W + P £ = ° > (8) 

9E 1 d . dib 

+ - — [rv r (E + P + p-f-v r + A nct = , 9 
at r Or Or 

AiP = AttGp , (10) 

where 

P = £ft, (11) 

i 

P = E^ = E^ T > ( 12 ) 

i i 

* = E (^r + ' <13> 

where p, n, f r , P, ^, E 1 , G, and k are the mass density, number density, radial velocity, gas pressure, 
gravitational potential, energy per unit volume, gravitational constant, and Boltzmann constant, 
respectively. The symbol A nct denotes the cooling function which represents the net energy loss 
rate per unit volume. The values with subscript i denote those of the i-th species. 

The number density of the i-th species, rij, is obtained by solving the following time-dependent 
rate equations, 

dx- 9 9 9 9 9 

= n H £ £ k jk XjX k + 4EEE hmnXlX m X n , (14) 
j=l k=l !=1 ra=ln=l 

where wh denotes the total number density of hydrogen nuclei and Xi = rii/nu is the relative 
number density of the i-th species. The reaction rate coefficients, kjk and ki mn , are given in Table 
1. 

In the models calculated in this paper, the gas temperature does not become far above 10 4 K. 
Thus the cooling is dominated by contributions from H at T ~ 10 4 K and H2 at T < 10 4 K. We 
therefore adopt the cooling rate in the form of A nct = Ah + Ah 2 +A c h em ) where Ah, Ah 2 , and A c h em 
denote contributions from H, H2, and chemical reactions, respectively. Ah includes the cooling of 
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Table 1. REACTION RATE COEFFICIENTS 



Reactions 



Rate Coefficients 



Reference 



(1) H + e^H+ + 2e 

(2) He + e -> He+ + 2e 

(3) He+ + e -» He++ + 2e 

(4) H+ + e -» H + hv 

(5) He+ + e -» He + /«/ 



5.85 x 10- n r - 5 exp(-157809.1/T)(l + T 5 

-1-11T-.-0.5 , 



0.7\-l 
0.5\-l 



(6) 


He ++ + e 


-> He+ + /«/ 


fc 6 r 




H + H -> 


H + H+ + e 


fc 7 - 


(8) 


H + e -> H" + 


/c 8 ; 


(9) 


H + H" - 


* H 2 + e 


fcg : 


(10) 


H + H+ -» H^ + /u/ 


fcio 


(11) 


H+ + H- 


* H 2 + H+ 


fell 


(12) 


H 2 + H- 


3H 


fcl2 


(13) 


H 2 + H+ - 


^H+ + H 


fcl3 


(14) 


H 2 + e^ 


H + H 


fci4 


(15) 


H 2 + e -» 


2H + e 


^15 


(16) 


H 2 + H 2 


-> 2H + H 2 


fcl6 


(17) 


H + e - 


-> H + 2e 


fel 7 


(18) 


H+H 


-> 2H + e 


fcl8 


(19) 


H - + H+ 


-» 2H 


feig 


(20) 


H" + H+ 


- H+ + e 


&20 


(21) 


H+ + e - 


-> 2H 


&21 


(22) 


H+ + H- 


-> H + H 2 


&22 


(23) 


3H -» H 2 


+ H 


&23 


(24) 


2H + H 2 


-> 2H 2 


&24 



-0.5 T3 -0.2 (1 + r6 0. 7) -l 



fcl 

fc 2 = 2.38 x 10- n T- a5 exp(-285335.4/T)(l +T 5 L 
fe 3 = 5.68 x 10- 12 T°- 5 exp(-631515.0/T)(l +T5 - 5 )- 1 
fc 4 = 8.40 x lO-^T- 1 / 2 ^- - 2 ^ + Tq - 7 )- 1 
A; 5 = 1.9 x 10- 3 T- 1 - 5 exp(-470000/T)(l +exp(-94000/T) 
+ 1.50 x i - 10 T-°- 6353 
3.36 x 10- 10 T" 
1.7 x 10- 4 fei 
1.0 x 10~ 18 T 

dex [-14.10 + 0.1175 logT 
-9.813 x 10- 3 ] (logT) 2 ] 
1.3 x 10~ 9 

dex [-8.78 + 0.113 logT 

-3.475 x 10- 2 (logT) 2 ] 
1.85 x 10- 23 T 18 
5.81 x lO" 16 

X (T /56200) M-6657 log(T/56200)] 

6.4 x 10" 



T < 1.5 x 10 4 K 

T > 1.5 x 10 4 K 
T < 10 4 K 

T > 10 4 K 

T < 6.7 x 10 3 K 



T > 6.7 x 10 3 K 



-10 



= 2.4 x 10~ 9 exp(-21200/T) 

= 2.7 x 10- 8 T- L5 exp(-43000/T) 

= 4.38 x 10- 10 T°- 35 exp(-102000/T) 

[see equation (5) in reference] 

= 4.0 x 10~ 12 exp(-43000/T) 

= 5.3 x 10~ 20 T 2 - 17 exp(-8750/T) 

= 7.0 x 10- 7 T-°- 5 
r 10 -8 T -0.4 

" \ 4 x 10^ 4 T- 1 - 4 exp(-15100/T) 

= 1.68 x 10- 8 (T/300)-°- 29 

= 5.0 x 10- 6 T-°- 5 

= 5.5 x lO^T- 1 

= k 23 /8 



T < 10 4 K 
T > 10 4 K 



1 
1 
1 
1 
1 

1 

2 



4 
5 
3 
6 
3 
5 
3 
3 
7 



7 
9 
9 



a The units of rate coefficients are taken to be cm 3 s 1 for two-body reactions and cm 6 s 1 for 
three-body reactions 

Note. - - (1) Cen 1992; (2) Palla et al. 1983; (3) Shapiro & Kang 1987; (4) Karpas, Anicich, 
& Huntress 1979; (5) Lepp k Shull 1983; (6) Hirasawa 1969; (7) Dalgarno & Lepp 1987; (8) 
Nakashima, Takayi & Nakamura 1987; (9) Palla et al. 1983 
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H atoms due to recombination, collisional ionization, and collisional excitation (Cen 1992; see also 
Black 1981). Ah 2 includes the cooling of H2 due to rotational and/or vibrational excitations (Lepp 
& Shull 1983; Haiman et al. 1996), collisional dissociation (Lepp & Shull 1983), and heating due 
to H2 formation (Shapiro & Kang 1987). A c h em includes the cooling due to chemical reactions in 
Table 1 (Shapiro Sz Kang 1987). For each contribution, we use the analytic formula expressed in 
each reference. [It is claimed that the cooling rate by Lepp & Shull (1983) is overestimated in 
relatively low densities (n# < lcm~ 3 ), compared to the rate by Hollenbach & McKee (1989). But, 
they are in a good agreement with each other in higher densities relevant to the present issue (Galli 
& Paha 1998).] 

When the H2 density exceeds ~ 10 10 cm~ 3 , the H2 line emission becomes optically thick 
(Palla et al. 1983). To take account of this effect, we modify the H2 line cooling with using 
the photon escape probability (Castor 1970; Goldreich & Kwan 1974): A^iine = /3AH 2 iine,thin and 
(3 = (1 — e~ TR ) /tr, where AH 2 iine,thin is the cooling function in optically thin regime, j3 is the photon 
escape probability, and tr is the Rossland mean opacity for the H2 line emission. To evaluate tr, 
we determined the level populations with the method of Palla et al. (1983). This type of (3 is 
applicable rigorously when the gas flows supersonically and its velocity monotonically changes in 
proportion to the radius (Castor 1970; Goldreich & Kwan 1974). However, the evolution has turned 
out to depend weakly on the form of (3. Using other types of (3, e.g., (3 = e~ TR , we have recalculated 
the evolution and reached basically the same results. 



2.2. A Model for a Filamentary Primordial Cloud 

In a gravitational instability scenario for the formation of the first structures, a cosmological 
density perturbation larger than the Jeans scale at the recombination epoch forms a flat pancake- 
like disk. This process has been extensively studied by many authors (e.g., Zel'dovich 1970; Sunyaev 
& Zel'dovich 1972; Cen & Ostriker 1992a, 1992b; Umemura 1993). Although the pancake formation 
is originally studied by Zel'dovich (1970) in the context of the adiabatic fluctuations in baryon or 
hot dark matter-dominated universes, recent numerical simulations have shown that such pancake 
structures also emerge in the CDM cosmology (e.g., Cen et al. 1994). Thus, the pancakes arc 
thought to be a ubiquitous feature in gravitational instability scenarios. 

In the bottom-up theory like a CDM model, the first collapsed objects at z ~ 10 — 100 are 
expected to have the masses of ~ 1O 5_7 M . In these clouds the gas is heated above T > 10 3-4 K 
by shock. Therefore, just after the shock formation, the cooling timescale is likely to be shorter 
than the free-fall one; the temperature of the pancake reduces to the value at which the cooling 
timescale is comparable to the free-fall one (e.g., Haiman et al. 1996; see also Yoneyama 1972). 
The temperature of the pancake is then estimated to be T ~ 200 — 1000 K for yn 2 = 10 -4 ~ 10 -3 . 
Thereafter, the pancake is likely to fragment into filamentary clouds which have nearly the same 
temperature as that of the parent pancake (Miyama et al. 1987a, 1987b; Uehara et al. 1996). In 
the following we describe a model of a filamentary cloud formed by fragmentation of the pancake. 
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As a model of a filamentary gas cloud, we assume an infinitely long cylindrical gas cloud for 
simplicity. At the initial state, we assume the gas to be quiescent and isothermal. We also assume 
that H~ , , He + , and He ++ do not exist at the start and the relative abundances of other species 
are spatially uniform for simplicity. The density distribution in the radial direction is set to be 



and 



Rm = J-^- , (16) 



where po and To are respectively the central density and the initial gas temperature, and / represents 
the degree of deviation from the equilibrium state. A model is specified by the following five 
parameters: po, Tq, f, the electron number fraction x e , and the H2 number fraction xn 2 . When 
/ = 1, the cloud is just in hydrostatic equilibrium; the density distribution accords with that of 
an isothermal filamentary gas cloud in equilibrium (Stodolkiewicz 1963). In this paper, we restrict 
the parameter / to / > 1 since we are interested in the evolution of the collapsing clouds. 

For the above model, the mass per unit length (line mass) is given by 

l = / 2irprdr = TrR m p c f 
Jo 

= *M° f = 2 .2 x 10 3 ( 1) A/opc- 1 ( -*°-) . (17) 

Note that the line mass of the equilibrium filamentary cloud depends only on the temperature. 
When the filamentary cloud forms through gravitational fragmentation of the pancake, its line 
mass is estimated to be I ~ 2p c [\ m H c i ~ 4kT/(pG), which is twice that of the equilibrium cloud; 
i.e., the typical value of / is evaluated to be / = 2 (Miyama et al. 1987a). Here p^, A m , and are 
the mass density of the parent pancake, the wavelength of the most unstable linear perturbation, 
and the half thickness of the pancake, respectively. We thus adopt / = 2 for a typical model. 



2.3. Model Parameters and Numerical Method 



As shown in § |2.2| , the initial state is specified by the parameters tiq{= po/p), 7b, /, x e , and 



xh 2 . We examine 60 models, by choosing the parameters no, To, and / to be no = 10 2 , 10 3 , 
10 6 cm -3 , T) = 10 2 , 5 x 10 2 , 10 3 , 10 4 K, and / = 2, 4, 10, respectively. For the parameters x e 
and sh 2) we adopt x e = 5 x 10 -5 and xn 2 = 10~ 4 for the models with To < 10 4 K. The value of 
x e = 5 x 10 -5 is adopted from the calculation of Peebles (1968) for the residual post-recombination 
ionization. [We have found that the results do not depend upon x e as far as x e ^ 10~ 7 , because the 
free electrons could quickly recombine to a level of x e < 10 -7 in the course of the collapse, (see also 
Haiman, Rees, & Loeb 1996, 1997).] The other abundances are determined by the conservation of 
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mass and charge. For the models with To = 10 4 K, we determine the abundances from the statistical 
equilibrium with e, H, H + , He, He + , and He ++ . 

The hydrodynamic equations (0) - (|l3|) are solved numerically by using the second-order upwind 
scheme based on Nobuta & Hanawa's (1998) method. This scheme is an extension of Roe's (1981) 
method to the gas having non-constant 7. [See Nobuta & Hanawa (1998) for more details and the 
test of the code.] The rate equations (|14|) are solved numerically with the LSODAR (Livemore 
Solver for Ordinary Differential equations with Automatic method switching for stiff and non-stiff 
problems) coded by L. Petzold and A. Hindmarsh. 

As for the boundary condition, we take the fixed boundary at r = R max , where -R max is the 
maximum radial coordinate in the computational domain. In all the models we have taken -R max to 
be about ten times greater than the effective radius of the filamentary cloud, f l / 2 R^\. The effect of 
the fixed boundary is very small. This is because the density is much lower near the outer boundary 
than at the center, i.e., p ^ I0~ 4 p c for all runs. (In fact, we have calculated the case with larger 
-Rmax and have confirmed that the numerical results are not changed.) 

The numerical grids are non-uniformly distributed so as to enhance the spatial resolution near 
the center. The grid spacing increases by 5 % for each grid zone with increasing distance from 
the center. As shown in §||, the characteristic scale shortens in the central high density region as 
the collapse proceeds. The spatial resolution thus becomes poor near the center. To compute the 
further contraction with the sufficient spatial resolution, we pursue the subsequent evolution with 
refining grids. Whenever the number of grid points becomes less than ~ 20 within the radius of 
the half-maximum density (p = 0.5p c ), we increase the number of grid points and then reposition 
them on the refined grids in the whole computation region (0 < r < -R max )- The physical variables 
at the new grid points are determined by linear interpolation. This technique allows us to pursue 
the dynamical evolution over more than the tenth order of magnitude in the central density. 



3. Numerical Results 

3.1. Clouds with Intermediate and High Initial Temperatures (To ^ 500 K) 

We have examined the evolution of various models with relatively high initial temperature 
(To ^ 500 K). As a result, the way of the evolution has turned out to be almost insensitive to the 
model parameters. 

As a typical example, we first show the evolution of the model with (no, To, /) = (10 4 cm -3 , 10 3 K, 2). 
Figure [l] shows the evolution of the temperature (T c ) and the mass fraction of H2 at the center 
(t/h 2 = 2xh 2 ) as a function of the central density n c . During the early evolution, H2 forms mainly 
through the H _ process [eqs. (|l|) and @]. Because of the effective cooling by the rotational tran- 
sitions, the temperature first descends promptly to T ~ 300 K. When the central density increases 
to n c ~ 10 6 cm -3 , the cloud collapses nearly isothermally, staying the temperature at T ~ 800 
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K until the central density reaches n c ~ 10 10-11 cm -3 . The temperature evolution can be then 
approximated by T c oc n^ 15 as far as n c < 10 8 cm -3 . When the density exceeds ~ 10 8 cm -3 , the 
hydrogen molecules form acceleratively through the three-body reactions [eqs. (||) and (||)]. After 
the central density increases to n c ~ 10 11 cm -3 , the hydrogen gas is almost completely processed 
into molecules. Then the cloud becomes optically thick against the H2 lines and the temperature 
rises up to ~ 1500 K. At this stage, the cooling time is 50 — 100 times longer than the free-fall 
time. Thereafter the collapse decelerates near the center and the cloud weakly oscillates around its 
quasi-static equilibrium state. 

In Figures ||a - ||d, the time variations of the distributions of the density, the temperature, the 
H2 mass fraction, and the contraction time are plotted at several dynamical stages. The contraction 
time is a dimensionless one, which is defined as t C ont( r ) = r l K-( r )%]) where tg is the free-fall time 
at the center, tg = (47rG/? c ) -1 / 2 . Since the gravitational force is perceptibly greater than the 
pressure force at the initial state, the cloud collapses nearly in a free-fall time. When the central 
density reaches n c ~ 10 8 cm ~ 3 (t = 1.434 x 10 6 yr), the shock arises around r ~ 3 x 10 2 pc, which is 
characterized by a jump in density and temperature. The shock surface separates the cloud into two 
parts which correspond to mutually different dynamical states, that is, a denser spindle and a diffuse 
envelope. During the contraction, the temperature in the spindle keeps nearly constant around 800 
K until the cloud becomes optically thick against the H2 lines. In the isothermal contraction phase, 
the outer parts in the spindle exhibit a power-law density distribution, p oc r~ 2 . It indicates that 
the spindle contracts in a self-similar manner. Actually, we have found a similarity solution, which 
is presented in the Appendix, and it has turned out that the newly found similarity solution well 
reproduces the numerical results. (See the Appendix for the further detail) At the shock front, the 
temperature rises up to T ~ 1500 K. The contraction time in the spindle is much shorter than 
that in the envelope. Thus the spindle collapses almost independently of the envelope. When the 
cloud becomes optically thick against the H2 lines at the center, the pressure force overwhelms the 
gravity, so that the second shock forms around r ~ 2 x 10 -4 pc. After this stage the contraction 
time of the spindle becomes much longer than % and the collapse promptly decelerates. 

Figure || shows the evolution of the line mass of the spindle as a function of the central density. 
We define the line mass as the density integrated over r up to the radius at which the density goes 
down to one tenth of the central density, 

rr(p=0.1p c ) 

l sp = / 2irrpdr . (18) 

Jo 

It is worth noting that during the contraction, the line mass of the spindle keeps a nearly constant 
value, which is ~ 1 x 1O 3 M pc -1 . This value is close to the line mass of an isothermal filamentary 
cloud in equilibrium for the temperature of 800 K, i.e., l = 2kT/pG = 0.9 x 10 3 M Q p C - 1 (T/800K). 
It implies that the gravitational force nearly balances with the pressure force in the spindle. 

We have found that the evolution is almost independent of the initial parameters, no, Tq, and 
/, for the models with Tq ^ 500 K. 
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3.2. Clouds with Low Initial Temperature (T ~ 100 K) 

In this subsection, we examine the evolution of models with relatively low initial temperature; 
To = 100 K. We find that the evolution in the case of low initial temperature depends mainly 
on the value of /. In Figures |]a and Qb, two models with (no, To,/) = (10 4 cm~ 3 , 10 2 K, 2) and 
(10 4 cm -3 , 10 2 K, 10) are compared. For both models the early evolution is similar; the clouds 
collapse adiabatically, T oc p 2//3 , because the H 2 cooling is not effective for T <^ 300 K. However the 
later evolution is quite different. For the model with a small /, the collapse almost ceases until the 
central density reaches 10 5 cnr 3 . In this case, the contraction time, which is nearly equal to the 
cooling time, is a hundred times as long as the free-fall time. Hence, the cloud becomes close to 
the hydrostatic equilibrium with no efficiency of cooling. On the other hand, for the model with a 
large /, the evolution is similar to that of the model with To ^ 500 K. When the cloud temperature 
exceeds 500 K, the collapse proceeds nearly isothermally due to the H 2 line cooling. In all the 
models with To = 100 K, it has turned out that the evolution proceeds in a similar way dependent 
upon the value of /. 



4. Fragmentation of Primordial Gas Clouds 

As shown in the previous section, if To ;> 500 K or / ;> 10, the filaments collapse quasi- 
statically, whereas the filaments would not continue to contract if To ^ 100 K and / <^ 2. In 
the latter case, the filament would be weakly unstable because the length of the realistic filament 
is shorter than the wavelength of the most unstable mode. Consequently, the filament is likely 
to contract along the axis without fragmentation. Also, as mentioned in § p.2j , the temperature 
of the initial filaments is estimated to be T ~ 200 — 1000 K in a CDM model. Hence, such low 
temperature clouds are hardly expected in a CDM scenario. Thus, here we focus on the former 
case, i.e., a collapsing filament, and consider the hierarchical fragmentation. 

The density of a collapsing filament is enhanced by more than eighth order of magnitude, with 
the constant temperature of T ~ 800 K, until the cloud becomes opaque. Thus, the filament could 
fragment into denser lumps. To diagnose the fragmentation process, we apply a linear theory to 
the numerical results and thereby estimate the length and mass scales of fragments. 

According to the linear theory of an isothermal hydrostatic filament, the perturbation grows 
most rapidly when its wavelength is about four times longer than the effective cloud diameter. 
The wavelength and growth rate of the most unstable perturbation depend on the density and 
temperature of the cloud. Therefore, they change according as the collapse proceeds. In our model 



of §3.1, the contraction timescale (t C ont) is longer by a factor of more than five than the free-fall time 
at the center (tg) after the central density reaches ~ 10 6 cm -3 . According to Inutsuka & Miyama 
(1992), when i cont ^ lOtg, the dispersion relation for the collapsing filament is approximated by 
that of the hydrostatic filament whose scale height is the same as the temporal one of the collapsing 
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cloud, and thus the growth of the density perturbation is approximated by 

5p(t)/p(t) = Aexp (ik z z - J iu(k z , t')d£\ , (19) 
during its linear growth phase, where k z and u denote the wavenumber and the frequency of 



the perturbation, respectively. We thus apply the above equation to our model of §3.1. For the 
growth rate (— iu), we use a fitting formula obtained by Nakamura, Hanawa, & Nakano (1993). 
A[= A(k z )] denotes the initial amplitude of the perturbation. We take A in the form of the power 
law A = A (k z /k z fi) p . When the amplitude of the unstable density perturbation becomes greater 
than unity (Sp/p ;> 1), the evolution of the perturbation becomes nonlinear and the cloud will 
breaks into pieces. It is expected that the evolution of the perturbation depends on the index p 
and the amplitude Aq. 

Assuming p = 1 or p = — 1, we have calculated the time evolution of the density contrast in the 
filamentary cloud. The results are shown in Figure ||. The solid curves show the density contrast 
at (1) t = yr (n c ~ 10 4 cm" 3 ), (2) 1.384 x 10 5 yr (10 7 cm" 3 ), (3) 1.477 x 10 5 yr (10 9 cm" 3 ), (4) 
1.500 x 10 5 yr (10 11 cm -3 ), and (5) 1.503 x 10 5 yr (3 x 10 12 cm -3 ), respectively. Here, we have taken 
(A ,k z>0 ) = (1.0 x 10 -3 , 1.26 x 10 2 pc -1 ). When n £ 10 10-11 cm -3 , the wavelength of the fastest 
growing perturbation is about a hundred times as long as the effective cloud diameter, because of 
the cumulated growth rate. Hence, the cloud is likely to tear into long filaments before the central 
density reaches n c ~ 10 11 cm -3 , if the fluctuations have enough initial amplitudes to enter nonlinear 
stages. The long filaments shrink further to form thinner filaments. As the shrink proceeds, the 
wavelength of the fastest growing perturbation shortens. Thus, the filaments may again tear into 
shorter and denser filaments. Consequently, hierarchical filamentary structures would form, if the 
initial amplitudes are larger than 0.001 as expected in a CDM cosmology. 

When the central density exceeds n c ~ 10 11 cm -3 , the cloud becomes opaque against the H2 
lines and the collapse quickly decelerates. Hence the hierarchical fragmentation would be termi- 
nated. Then the wavelength of the fastest growing perturbation becomes A, ~ 2 x 10 -3 pc, which 
is comparable to that expected from the linear theory of an isothermal filament in equilibrium, 
as shown by the growth rate in the ending stage in Figure 5. Then, the cloud is most unstable 
against the perturbation whose wavelength is about four times as large as the effective cloud di- 
ameter. This wavelength is nearly independent of the power index p as shown in Figure & As 
a result, the filament is likely to finally fragment into dense cores, the typical mass of which is 
A*J S p ~ (2 x 10 -3 pc) x (1.5 x 10 3 Mqpc -1 ) ~ 3Af Q . This is the lowest mass of fragments, which pro- 
vides the minimum core mass for the first-generation stars. The mass possibly increases by accreting 
the ambient gas. If all the ambient gas within one wavelength of the fastest growing perturbation 
accretes onto the core, then the core mass increases to A*/o ~ 16(//2)(To/1000 K)M@. Thus, the 
mass of the first stars is expected to be in the range of 3M© M <^ 16(//2)(To/1000 K)M & . It 
should be noted that the lower mass as well as the upper mass is not sensitive to the amplitude 
of the power-law spectrum, because the lower mass is basically determined by the micro process of 
the H2 cooling and the upper mass depends only upon the initial temperature. 
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5. Conclusions 

In a wide range of the parameter space, we have numerically explored the thermal and dynam- 
ical evolution of filamentary primordial gas clouds with including nonequilibrium process for the 
H2 formation. We have found that for the models with To ^ 500 Kor/^ 10, the collapse proceeds 
nearly isothermally due to the H2 line cooling, staying the temperature at T ~ 800 K. During the 
contraction, the cloud structure is divided into two parts, i.e., a spindle and an envelope. The line 
mass of the spindle keeps a nearly constant value of ~ 1 x 1O 3 M0 pc _1 during the contraction from 
10 6 cm -3 to 10 11 cm~ 3 . The outer part in the spindle exhibits a power-law density distribution 
as p oc r -2 . This behavior is well reproduced by a newly found self similar solution. When the 
central density reaches n c ^ 10 11 cm -3 , the cloud becomes opaque against the H2 line cooling and 
reaches a quasi-static equilibrium state with keeping the line mass. On the other hand, for the 
models with To ~ 100 K and / <J 2, the cloud contraction is much slower because the H2 cooling is 
not effective. The contraction immediately ceases. In a CDM scenario, the former case (collapsing 
filaments) would be more probable. 

Applying a linear theory for the gravitational instability of collapsing filaments, we find that the 
spindle is unstable against fragmentation during the quasi-static contraction. The wavelength of the 
fastest growing perturbation (A m ) lessens as the collapse proceeds. Thus successive fragmentation 
could occur. When the central density exceeds n c ~ 10 11 cm -3 , the successive fragmentation 
may be suppressed since the cloud becomes opaque against the H2 lines and the collapse promptly 
decelerates. Resultingly, the minimum value of A m is estimated to be ~ 2 x 10 -3 pc. The typical 
mass of the first stars is then expected to be ~ 3M , which may grow up to 16Mq by accreting the 
diffuse envelope. The present results may be relevant to the early evolution of primordial galaxies 
and the metal enrichment of the intergalactic space. 

We are grateful to Drs. T. Nakamoto, and H. Susa for stimulating discussion. We also thank 
an anonymous referee for valuable comments which improved the paper. This work was carried out 
at the Center for Computational Physics of University of Tsukuba. This work was supported in 
part by the Grants-in Aid of the Ministry of Education, Science, and Culture (09874055, 09740171, 
10147205). 
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A. Similarity Solutions of A Polytropic Filamentary Cloud 

As mentioned in §3, the spindle seems to collapse self-similarly during the nearly isothermal 
contraction phase from n c ~ 10 6 cm -3 to 10 11 cm -3 ; e.g., the density distribution of the spindle is 
characterized by the power-law distribution of p oc r~ 2 . But this behavior is different from that of 
the similarity solution (p oc r -4 ) of an isothermal filamentary cloud that is known so far (Miyama 
et al. 1987a; see also eqs. [ A18 | and [ A19f )■ In this appendix we seek another type of similarity 



solutions in more generic equations of state. 

We consider an infinitely long cylindrical cloud in which the density is uniform along the axis. 
We assume the gas to be polytropic; P = Kp^, where K is constant. The equation of motion and 
the continuity equation are then described as 

dv r dv r 2GI ldP _ 
dt r dr r p dr 

— + 2irrpv r = , A2 
at 

and 

dl 

— - 27rrp = , (A3) 
where / denotes the line mass (mass per unit length) contained within a radius r. 
We now look for a similarity solution of the form 



r 

at 



(A4) 



*,Q-£L, (A5) 

v r (r, t) = av r (x) , (A6) 

l(r,t) = ^!(x) , (A7) 

and 

a = t l -^^K{2TrGY^ , (A8) 

where the physical variables with a tilde denote those in the similarity coordinate, x. Substituting 

(A9) 



equations (|A4[) — (A8) into equations ([Al]) — (|A3|) , we obtain 

dp Y(x,p,v r ) 



dx X(x, p, v r 
dv r Z(x,p, v r 



ix X(., A *)' (A10) 

X(x,p,v r ) = [«,. - (2 - f)x] 2 - TP 1 '" 1 , (All) 
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Y(x, p, v r ) = - [v r - (2 - j)xY 



7-1 



+ (2x 



(7 - l)pv r 



(A12) 



and 



Z(x, p, v r 



(2 - 7 )x] 2 ^- + (7 - l)v~r[vr - (2 - 7)x] + ^ 



,7-1 



-iy r -2x). (A13) 



7 — 1 x 

Note that equations (|A9[ ) through ( A13| ) have a singularity at 7 = 1 (A polytropic sphere has a 
singularity at 7 = 4/3. This is fundamentally the same as the dynamical stability condition of the 
sphere.). 

One of the solutions of equations (K9) through ( A13| ) is expressed as 

v r = 



and 



(2-7) 2 



1 

7-2 



X 2 -T 



(A14) 



(A15) 



.27(1-7). 

In this solution the gas is at rest and the density diverges at the origin. This solution corresponds 
to a singular equilibrium solution for a polytropic filamentary cloud. 

A regular solution at x = has an asymptotic behavior, 

dhxp 2 



and 



dlnx 
d In v r 



2(7 



•7 
-1) 



(A16) 



(A17) 



d In x 2 — 7 

Thus, when the value of 7 is nearly equal to unity, the density and velocity distributions follow 
p — > r~ 2 and v r — ► const, in the outer region of |x| S> 1. This behavior is quite similar to that of a 
sphere (Larson 1969; Penston 1969). Note that this similarity solution exists only for 7 < 1. This is 
essentially the same as the dynamical stability condition of the polytropic cylinder. This similarity 
solution has different characteristics from a singular case of 7 = 1. In the similarity solution of 
7 = 1, the density and velocity distributions are given as (see Miyama et al. 1987a) 

4 



P 



(1 + x 



2\2 



and 



(A18) 
(A19) 



which are proportional to x 4 and x, respectively, in the region of \x\ ^> 1. 



In Figure ^ the similarity solutions of 7 = 0.9 and 0.99 are compared with that of 7 = 1. The 
thick and thin solid curves represent the density and velocity distributions of the similarity solutions 
with 7/I, respectively. The dashed curves are the similarity solution of 7 = 1 (Miyama et al. 



-15- 



1987a). The newly found similarity solution seems to well reproduce the numerical results during 
the nearly isothermal contraction phase; the density distribution in the outer part in the spindle 
is nearly proportional to r~ 2 and the velocity distribution is proportional to r near the center [see 
Fig. ||. Recently Kawachi & Hanawa (1998) also studied the gravitational collapse of a polytropic 
cylinder and confirmed that the evolution of the cylinder converges to the similarity solution.]. It 
should be stressed again that whether 7 is precisely unity or not transforms conclusively the self- 
similar manner. The present similarity solution for 7 ^ 1 is likely to be applicable for the realistic 
situations. 
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Fig. 1. — Evolution of the temperature (T c , solid curves) and the H2 abundance (yu 2tC , dashed 
curves) at the center for the model with (no, Tq,/) = (10 4 cm -3 , 10 3 K, 2). The first prompt 
decrease of the temperature is due to the H2 cooling. Afterwards, until the central density reaches 
~ 10 11 cm -3 , the temperature keeps a nearly constant value of ~ 800 K over the fifth order of 
magnitude in density. The H2 abundance steeply rises around n c ~ 10 9 cm 3 and reaches ?/h 2 ~ 1 
at the stage at which n c ~ 10 11 cm -3 . At n c ~ 10 11 cm -3 , the cloud becomes optically thick against 
the H2 lines, and consequently the temperature rises to ~ 1500 K. Thereafter the cloud contraction 
decelerates near the center and the cloud weakly oscillates around its quasi-static equilibrium state. 
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Fig. 2. — Time variations of the distributions for (a) the density, (b) the temperature, (c) the H2 
abundance, and (d) the contraction time are shown for the same model as Fig. 1, i.e., (no, Tq, f) = 



(10 4 cm 3 , 10 3 K, 2). They are plotted at the stages at which (1) t = yr (n c = 
t = 1.285 x 10 6 yr (n c = 10 6 cm" 3 ), (3) t = 1.434 x 10 6 yr (n c = 10 8 cm" 3 ), (4) t 
(n c = 10 10 cm" 3 ), and (5) t = 1.502 x 10 6 yr (n c 



10 4 ) cm" 3 , (2) 
= 1.496 x 10 6 yr 



10 12 cm" 3 ) 
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n c (cm 3 ) 

Fig. 3. — Evolution of the line mass of the spindle in units of M & pc -1 for the model with 
(no,?o, /) = (10 4 cm -3 , 10 3 K, 2). The line mass of the spindle keeps a constant value, which is 
~ 1O 3 M0 pc _1 until the cloud becomes optically thick against the H2 lines. This value accords 
with that of an equilibrium isothermal filamentary cloud with the temperature of T ~ 800 K. It 
indicates that in the spindle the gravitational force nearly balances with the pressure force. 
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Fig. 4. — Same as Figure |T| but for the models with (a) (no, To,/) = (10 4 cm~ 3 , 10 2 K, 2) and 
(b) (10 4 cm~ 3 , 10 2 K, 10). For both models the early evolution is similar; the clouds collapse 
adiabatically, T oc p 2 ^ 3 , because the H2 line cooling is not effective for such low temperature. For 
the model with a small / the collapse almost ceases until the central density reaches 10 5 cm~ 3 . On 
the other hand, for the model with a large /, the evolution is similar to that of the model with 
To ^ 500 K. When the cloud temperature exceeds 500 K, the collapse proceeds nearly isothermally 
due to the H2 line cooling. 
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Fig. 5. — Evolution of density perturbations for the model with (no, To, /) = (10 4 cm" 3 , 10 3 K, 2). 
The abscissa and ordinate denote the wavelength and amplitude of the density perturbation, 
respectively. The initial spectrum of the density perturbations is assumed to be in the form 
of A = Ao(k z /k z fi) p . The index p is (a) p = 1 or (b) p = —1. We take (Ao,k Zt o) = 
(1.0 x 10~ 3 ,1.26 x 10 2 pc _1 ). The solid curves represent the amplitude of the density pertur- 
bations at (1) t = yr (n c ~ 10 4 cm" 3 ), (2) 1.384 x 10 6 yr (10 7 cm" 3 ), (3) 1.477 x 10 6 yr (10 9 
cm" 3 ), (4) 1.500 x 10 6 yr (10 11 cm" 3 ), and (5) 1.503 x 10 6 yr (3 x 10 12 cm" 3 ), respectively. For 
comparison, we show the square of the growth rate at the stage (5) with the dashed curves. 
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Fig. 6. — Similarity solutions of filamentary clouds with 7 = 0.9 and 0.99. The thick and thin solid 
curves show the density and velocity distributions in the similarity coordinates, respectively. For 
comparison, we denote the density and velocity distributions of the isothermal similarity solution 
with the thick and thin dashed curves, respectively. For the similarity solutions with 7 = 0.9 and 
0.99, the density distribution is nearly proportional to r~ 2 and the velocity converges to a constant 
value in the outer part. 



